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ABSTRACT 

The stable clustering hypothesis is a key analytical anchor on the nonlinear dy- 
namics of gravitational clustering in cosmology. It states that on sufficiently small 
scales the mean pair velocity approaches zero, or equivalently, that the mean number 
of neighbours of a particle remains constant in time at a given physical separation. N- 
body simulations have only recently achieved sufficient resolution to probe the regime 
of correlation function amplitudes ^ ^ 100 — 10^ in which stable clustering might be 
valid. In this paper we use N-body simulations of scale free spectra P{k) oc A:" with 
— 2 < n < and of the CDM spectrum to apply two tests for stable clustering: the time 
evolution and shape of ^(x, t), and the mean pair velocity on small scales. We solve the 
pair conservation equation to measure the mean pair velocity, as it provides a more 
accurate estimate from the simulation data. For all spectra the results are consistent 
with the stable clustering predictions on the smallest scales probed, x < 0.07 Xni{t), 
where Xni(t) is the correlation length. The measured stable clustering regime corre- 
sponds to a typical range of 200 ^ £, ^ 2000, though spectra with more small scale 
power (n ~ 0) approach the stable clustering asymptote at larger values of ^. 

We test the amplitude of £, predicted by the analytical model of Sheth & Jain 
(1996), and find agreement to within 20% in the stable clustering regime for nearly all 
spectra. For the CDM spectrum the nonlinear ^ is accurately approximated by this 
model with n ~ —2 on physical scales ^ 100 — 300/i~^kpc for cjg = 0.5 — 1, and on 
smaller scales at earlier times. The growth of ^ for CDM-like models is discussed in the 
context of a power law parameterization often used to describe galaxy clustering at 
high redshifts. The growth parameter e is computed as a function of time and length 
scale, and found to be larger than 1 in the moderately nonlinear regime - thus the 
growth of ^ is much faster on scales of interest than is commonly assumed. 
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1 INTRODUCTION 

The stable clustering hypothesis is one of the few analytical 
handles on the deeply nonlinear regime of gravitational clus- 
tering. It states that the mean relative velocity of particle 
pairs in physical coordinates is zero. Hence in comoving co- 
ordinates the mean relative velocity exactly cancels the Hub- 
ble recession velocity between pairs of particles. As we will 
see below, this is equivalent to the statement that the mean 
number of neighbours of a particle in physical coordinates 
remains constant in time. The stable clustering hypothesis 
invokes the physical picture of a virialized cluster which has 
separated out from the expanding background, and is nei- 
ther expanding nor contracting. Since on small scales, any 
statistical measure is dominated by the contribution from 
dense clusters, the clustering should statistically be stable. 



This hypothesis leads to predictions for the evolution in 
time of the autocorrelation function ^. In the case of scale 
free spectra which display self-similar scaling, the slope of 
^{x) in the small scale regime is predicted as well. For ^ ^ 1 
it is extremely difficult to make analytical approximations 
to the growth of clustering. Hence the stable clustering hy- 
pothesis has been very useful in relating the shape of ^ and 
the power spectrum P{k) in the nonlinear regime to the 
initial spectrum. Peebles (1974) and Gott & Rees (1975) 
implemented the first such applications. The widely used 
property of hierarchical scaling of the higher order moments 
of the density in terms of its second moment also derives 
from the dynamics of stable clustering. As emphasized in 
the next section, this property requires additional assump- 
tions of stability at each order in the distribution. 

N-body simulations are the ideal tool to test stable clus- 
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tering. Efstathiou et al. (1988) carried out a pioneering study 
of self-similar evolution, and tested the small scale £, and v 
for stable clustering. Their data were consistent with the sta- 
ble clustering slope for ^ on the smallest scales, but their 32^ 
particle simulations lacked the resolution to give any defini- 
tive conclusion. Recently Padmanabhan et al. (1995) and 
Colombi, Bouchet & Hernquist (1995; hereafter CBH), have 
tested the stable clustering predictions for ^,v, and in the 
latter case for the hierarchical form of the higher moments 
5*3, Si and S5 as well. While Padmanabhan et al. (1995) 
found departures from stable clustering in their n = 1 sim- 
ulations, CBH verified the stable clustering prediction in 
their tests of ^. However, they found small departures from 
the hierarchical relation for the higher moments. 

In this paper we use P^M (particle-particle/particle- 
mesh) simulations with 100^ — 144^ particles to test stable 
clustering for power law spectra with —2 < n < 0, and 
for the CDM spectrum. We measure ^ and v, for which we 
use the pair conservation equation to obtain estimates with 
greater accuracy than allowed by direct measurements. Our 
approach and numerical resolution is similar to that of CBH, 
with the advantage that we have 4 — 10 times as many par- 
ticles. We do not however test for the higher moments as 
they do. While stable clustering predicts the slope of its 
amplitude can be approximated by additional assumptions 
as done in the analytical model of Sheth & Jain (1996). We 
test their predictions against the measured amplitude in the 
N-body data. We also test for stable clustering in two high 
resolution CDM simulations and identify the range of scales 
and epochs on which it is an adequate approximation. 

We begin in the next section with the relevant BBGKY 
equations through which the stable clustering hypothesis 
leads to the predictions for ^ and the higher moments. In 
Section 3 we outline a secondary infall model which can 
be connected with the nonlinear form of ^, and discuss the 
dynamical effects which might invalidate stable clustering. 
Section 4 contains a description of the N-body simulations, 
and an analysis of the effects of limited numerical resolution 
on small scales. Section 5 presents the main results for ^, 
with the mean pair velocity results in Section 5.1, and re- 
sults for the CDM spectrum in Section 5.2. We conclude in 
Section 6. 



2 THE BBGKY HIERARCHY EQUATIONS 

The stable clustering hypothesis yields the form of the non- 
linear correlation function when applied to the moments of 
the BBGKY hierarchy equations. The presentation below 
follows that of Peebles (1980, Sections 71 and 73). 

The second equation of the BBGKY hierarchy of equa- 
tions, when integrated over momenta, yields the pair conser- 
vation equation. This equation connects the rate of change 
of the autocorrelation function ^{x,t) to the mean pair ve- 
locity in comoving coordinates, v(x,t) = ax: 



0. 



(1) 



Integrating the above equation over a sphere of radius x 
gives 



d_ 

dt 



dx'x''^(l+i) 



where n{t) cx a~'^ is the number density in physi- 
cal coordinates. From the definition of ^ it follows that 
n{t)a'^ATTx'^dx\l'\-^{x,t)] is the typical number of neighbours 
in a spherical shell of thickness dx at a distance x from a 
particle. Thus the first term in equation ^ is the rate of 
change of neighbours inside a sphere of radius x, while the 
second term is the flux of particles through the surface of 
the sphere - hence the name pair conservation equation. 

The stable clustering hypothesis is the statement that 
the relative pair velocity in comoving coordinates, ax, ex- 
actly cancels the Hubble velocity, ax, so that the mean pair 
velocity in physical coordinates d{ax)/dt = 0. On substitut- 
ing this assumption. 



v{x, t) 



-ax. 



(3) 



in (|l]), the functional form of ^(k, i) which satisfies the equa- 
tion becomes restricted to 



I + e. = a\f [a{t)x\ , 



(4) 



and we can approximate 1 ^ ~ ^, as the stable clustering 
hypothesis is applied on very small scales were ^ is at least 
of order 100. Equation (^) thus fixes the growth of ^ in 
physical coordinates: ^{r,t) oc a"^. Physically it means that 
the typical number of neighbors per unit volume, at small 
enough separations r from a particle, remains constant in 
time, as it is ~ n{t)^{r,t), and n{t) oc a"'' as noted above. 

If the universe is spatially flat, with Q = ilmatter ~ 1, 
and the initial spectrum of perturbations is a pure power law 
P{k) (X fc", the subsequent evolution is expected to be self- 
similar (Peebles 1980, Section 73). Then ^{x,t) and P{k,t) 
take the special functional form; 



^{x,t) =i{x/xoa°') ; P{k,t) ^ a^^k^^ P{ka" /ko) 



(5) 



where a — 2/(3 -I- n); ko, xq are constants which must be 
determined from the initial conditions; and P, ^ are unspec- 
ifled dimensionless functions. It is easy to verify that the 
linear spectrum Pi{k,t) oc a^k" is consistent with the func- 
tional form of equation (^). The linear mean squared den- 
sity contrast on length scale k~^ scales as a'^k~^''^'^"\ hence 
equating it to unity gives the same scaling, k oc a~'^^^^^"^ 
as contained in equation (^). This unique scaling in time 
of characteristic length scales is the essence of self-similar 
evolution. 

Equating the functional forms of ^ given by equations 
(^) and ^ restricts the t and x dependence of ^ to be a 
pure power law. The solution for ^, and similarly for P, is: 



C(x,t) oc ix/a°')-^ ■ P{k,t) oc a(J+^fcT5+^ 



(6) 



+ 4™ X n(l + C) »^ = , (2) 



where 7 = (9 -f 3n) /(5 4- n). Thus the stable cluster- 
ing assumption, in conjunction with the requirement of 
self-similarity, fixes both the asymptotic growth and the 
shape of the power spectrum/correlation function. While 
self-similarity is an idealization as it requires the initial spec- 
trum to be a pure power law, it can be a useful approxima- 
tion to a more realistic description, and as we shall see it 
is very useful for testing the dynamical validity of stable 
clustering. 

For completeness we re-write the pair conservation 
equation ^ in terms of the mean interior correlation func- 
tion, which is the integral of the correlation function over a 
sphere, divided by the volume of the sphere: 
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ax,t) = — / dxx^ax,t). (7) 

2^ Jo 

Using this definition and re-arranging terms in equation 
gives 



— V 



3(1 + da ax ~ 

where Hr is the Hubble velocity on the physical scale r = ax. 
In Section 4 we shall use this form of the pair conservation 
equation to measure the mean pair velocity v from the N- 
body simulations. 

2.1 Hierarchical relation of higher moments 

The hierarchical form of the 3-point function ("123 = 
C('"i2, r23, rai) and of the higher moments can also be related 
to the stability assumption. Here we outline the argument 
which motivates the hierarchical form using self- similarity 
and stable clustering. The purpose of this exercise is to point 
out that additional assumptions of stability at every order 
in the particle distribution are needed to obtain the hierar- 
chical form at that order. These go beyond the assumption 
—v/Hr — 1 used to obtain the stable clustering shape of ^. 

Following Peebles (1980), consider the equation of con- 
servation of triplets of particles, obtained by integrating 
over momenta the 3rd equation of the BBGKY hierarchy. 
In terms of the relative velocities of particle pairs, it can be 
written as 

Oh 

+ ■ ^la'ia + ■ wsih-i = , (9) 

at 0x12 oxii 

where /la = 1 + C12 + C23 + Cai + C123 , and Wij is the mean 
pair velocity of particles i,j given the position of the third 
particle in the triangle formed by the triplet. This definition 
makes W12 a different statistic from V12: the former is a 3- 
point statistic as it is obtained by integrating over momenta 
the 3-point distribution function, while the latter is a 2-point 
statistic. 

In equation if we assume that 1012 — —ax\2 and 
= — aa;3i on small scales, then we are led to the func- 
tional form h'i = a^h{axi2,ax23). Combined with the con- 
straint of self-similarity this leads to a power law solution 
for hs and therefore for the three point function: oc x~^^ , 
where x is the size of the triangle, and 7 is the stable clus- 
tering slope of ^ defined following equation (^. It is thus 
consistent with the form 



C123 = Q (C12C23 + C23C31 + C31C12) 



(10) 



which is the widely used hierarchical form for the 3-point 
function. This argument can be generalized to the 4- and 
5-point functions, and has been found to agree with results 
of N-body simulations as well as observations of galaxies. 

The argument leading to the hierarchical form for C, 
makes an assumption about the 3-point particle distribu- 
tion. Likewise for the A'^-th order correlation function, the 
assumption required is that the mean relative velocity of par- 
ticle pairs, given the presence of A'^ — 2 other particles, be 
zero. Therefore additional conditions of stability are needed 
at each order of the hierarchy to get the analogues of equa- 
tion (|l^) for higher order moments. These conditions clearly 
go beyond the assumption v = —ax which led to equation 
(a) for the second moment ^. The tests performed in this 



paper relate to 5, and therefore do not necessarily verify the 
assumptions leading to the hierarchical form of higher mo- 
ments. This is also relevant in comparing our results with 
those of CBH, who appear to find small departures from the 
hierarchical relation for the S3 , S4 and S5 parameters (while 
their results for ^ are consistent with ours as discussed in 
Section 5). 



3 ANALYTICAL MODELS 

For convenience in modelling, small scale clustering can be 
thought of in two somewhat distinct ways. One is the hi- 
erarchical growth of structure by the continuous mergers 
of smaller clumps. N-body simulations have highlighted the 
importance of mergers, but the detailed dynamics is very 
difficult to model analytically. The other model, known as 
secondary infall, visualizes a more gradual and spherically 
symmetric accretion of matter on to initial density peaks 
leading to the formation of a halo. Simulations show that 
spherically symmetric accretion rarely occurs in the forma- 
tion of typical halos. But one might imagine that the average 
properties of halos are well approximated by the outcome of 
such a secondary infall model. 

The advantage of the secondary infall picture is that it 
is analytically tractable, and the density profile of a halo re- 
sulting from a given initial profile can be calculated. Gunn 
& Gott (1972) and Gott (1975) made the first calculations 
of spherical accretion models, while Fillmore & Goldreich 
(1984) and Bertschinger (1985) obtained detailed solutions 
for the profile arising from the late time behavior of par- 
ticle orbits. Hoffman & Shaham's (1987) inffuential paper 
applied these results to the structure of halos formed in a 
cosmological setting. The secondary infall solution for the 
density profile can also be connected to the nonlinear form 
of 5 (Padmanabhan et al. 1995; Padmanabhan 1996; Sheth 
& Jain 1996). Here we outline this connection, and comment 
on its implications for the validity of stable clustering. This 
section is a bit of a detour from the main body of the paper, 
and some readers may wish to skip to sub-section 3.2 which 
provides a summary of the section. 

Consider a spherically symmetric distribution of coUi- 
sionless particles with an initial density profile 

^{X,ti) ocx^'^ . (11) 

p 

Note that k here is 3e in the notation of Fillmore and Goldre- 
ich (1984) as they used the mass as the variable on the RHS. 
Let these particles be assigned velocities such that they are 
in Hubble fiow at some initial time ti. The trajectories of 
these particles will follow radial orbits, and since they are 
in overdense regions, they will eventually stop expanding 
and collapse towards the center of mass. The point at which 
particles with an initial radius Xi first stop expanding is 
known as turnaround, and the turnaround radius denoted 
xta{xi) depends on Xi. After turnaround, since the particles 
are taken to be coUisionless, they execute oscillations about 
the center of mass. 

Fillmore & Goldreich obtained self-similar solutions for 
the final distribution of these particles by using the adia- 
batic invariance of the mass M{x, t) within a given radius x. 
Consider the orbit of a particle with maximum radius much 
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smaller than the current turnaround radius. By assuming 
M{x,t) to vary sufficiently slowly that it could be treated 
as a constant over the timescale of oscillation of this particle, 
they solved for the variation in time of the maximum radius 
of the particle in terms of its turnaround radius. This in turn 
allowed them to obtain the asymptotic (late time) density 
profile in terms of the initial density profile. The result is: 



{X,tasym) OC X <! + «) 

P 



K>2. 



(12) 



The shape of the final density profile in the secondary 
infall model is qualitatively different if k < 2, i.e. if the 
initial density profile is sufficiently shallow. The dynamical 
origin of this difference is that for k < 2, at late times the 
mass within a radius x is dominated by the contribution 
from particles with maximum radii > x. In contrast, for 
K > 2, the mass within x is dominated by particles with 
maximum radii < x. The consequence of this difference is 
that for K < 2, particle orbits on small scales keep shrinking 
due to the increasing mass contributed by large scales. The 
late time density profile approaches a form that does not 
depend on the initial profile. The asymptotic density profile 
is given by 



^Ax t 

p 



) OC X 



K.<2. 



(13) 



The above discussion relies on purely radial orbits. The 
sharp transition at k = 2 does not occur if one allows the 
particles to have some angular momentum. White and Zarit- 
sky (1992) have studied the density profiles arising from an 
initial distribution of angular momentum in which the eccen- 
tricity is constant as a function of radius. In such a case, the 
transition to a regime of density profiles does not occur 
at all, and the form of equation (112) is always valid. More 
generally one might find that a finite angular momentum 
causes the density profiles to lie somewhere in between the 
shapes of equation ( p^ and (^^ for k < 2, and to maintain 
the shape of equation (^2|) for k > 2. 



3.1 From secondary infall to stable clustering 

We are now in a position to apply the secondary infall model 
to the autocorrelation function ^ in an Einstein-de Sitter cos- 
mology. The hope is that what works for an isolated halo also 
applies statistically to the full matter distribution with an 
appropriate matching of initial profiles. Consider the iuCTe- 
dients that led to the stable clustering profile of equation dq). 

(i) The pair conservation equation, which is simply a kine- 
matical expression of the conservation of particle number. 

(ii) Self-similarity, namely that characteristic length scales 
grow as a power law of time: x oc a{t)^^''^'^"\ (iii) Expand- 
ing coordinates, implicit in the use of the comoving spatial 
coordiate x. 

Each of these ingredients is present in the secondary in- 
fall model described above. Thus stable clustering relies on 
such a limited aspect of the dynamics of the full distribution, 
that the model of an isolated spherically symmetric object 
can satisfy the conditions under which it is applied. To con- 
nect the two models, we need to relate the initial density 
profile of the secondary infall model to the initial spectrum 
of the cosmological distribution. Unfortunately there is no 
rigorous way to make this connection. We show below that 



the stable clustering slope for ^ is recovered if the initial 
density profile is taken proportional to the rms smoothed 
density fiuctuation of the cosmological distribution (Pad- 
manabhan et al. 1995) which is given by 

^(:r,tO OC x-'^ OC (14) 
P 

The Fillmore-Goldreich result for the final profile given 
by equations (jl^) and ( p^ ) can now be directly adapted, 
with the substitution « = (3 + n)/2. Thus the transition 
value of K = 2 corresponds to n = 1, and the asymptotic 
shape of ^ is 



(9 + 3") 

^{X,tasym) OC X <5 + ") 



n > 1 , 



(15) 



corresponding to equation (^), and, 

^{X,tasym) OC x"^ ; Tl < 1 , (16) 

corresponding to equation (p^. 

Equation (^^ shows one way to connect the stable clus- 
tering profile to the nonlinear density profile of an isolated 
spherical halo. Since the choice of the initial profile was made 
in hindsight, in order to recover the desired shape of ^, this is 
not intended to be a "derivation" of stable clustering. What 
it does in fact demonstrate is that there could be a regime in 
which stable clustering is invalid. With the particular initial 
profile chosen here, stable clustering holds only if the initial 
spectral index n > 1 - a range of very little interest for re- 
alistic spectra like the CDM spectrum, since for all scales 
of interest n ranges between 1 and —3 (or lower for spec- 
tra with some hot-dark matter). Taken at face value this 
result suggests that for spectra of interest stable clustering 
is invalid, and that the nonlinear ^ should take the universal 
shape ^ OC x~^, independent of the initial spectrum. 

However, as discussed above, the presence of nonradial 
orbits lowers the transition value of k and of n, and possibly 
completely eliminates the regime of an x~^ profile. The key 
parameter therefore becomes the degree of eccentricity of 
the orbits - if the eccentricity stays roughly constant, or 
increases with the size of the orbit, then infalling particles 
will not contribute to the mass inside a given x and orbits 
will not asymptotically shrink. In such a case, the profile 
will retain the stable clustering prediction of equation (p^. 

One can also consider a more detailed model in which 
pairs with small separation x are taken to belong to stable 
halos, so that £^{x) arises from the convolution of the den- 
sity profiles of such halos (Sheth & Jain 1996). The stable 
clustering shape of ^ can be obtained from the secondary in- 
fall result with an appropriate initial halo profile, but again 
there is no rigorous justification for this initial profile. 

3.2 Summary 

The lesson from the secondary infall models for ^ is that 
they show how the stability assumption can be invalidated in 
cases where the initial profile is sufficiently shallow. The key 
dynamical factor is the distribution of angular momentum. 
Particles with large orbital radii need to have sufficient an- 
gular momentum, else their contribution to the mass inside 
arbitrarily small radii causes small scale orbits to asymptot- 
ically shrink. This corresponds to a net infall pair velocity, 
and could lead to a universal shape for ^, independent 
of the initial n for spectra of interest. 
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In the formation of halos a more important process ap- 
pears to be the hierarchical merging of smaller halos. The 
influence of merging on small scale dynamics can depend 
critically on the size of halos which contribute most of the 
pairs at small separations. If most of the pairs with sepa- 
ration X belong to the cores of halos of radii much larger 
than a;, then ongoing mergers of roughly equal mass halos 
might cause the cores to expand and lead to a net outflow 
pair velocity. On the other hand if members of pairs come 
from halos with radii of order, or smaller than, x, then the 
merging of such halos could lead to net infall pair veloci- 
ties. These considerations highlight some mechanisms which 
could violate the stability assumption, and cause ^ to be- 
come either steeper or shallower than the stable clustering 
shape. But the limitations of analytical modelling leave N- 
body simulations as the only realistic means of testing for 
stable clustering. 



4 N-BODY SIMULATIONS 

This section provides some details of the N-body simulations 
used to obtain the results presented in the next section. Six 
different siirmlations arc analyzed in this paper. Two are for 
the standard CDM spectrum. The other four are for power- 
law spectra with n = 0, —1, —1.5 and —2 with Q, = 1. All 
the simulations were performed using high resolution P^M 
codes. Table 1 shows some important parameters for the 
different simulations: the total number of particles Npart, 
the Plummer force softening parameter e, the size of the 
PM mesh, and the smallest comoving scale Xmin on which ^ 
can be reliably measured. 

For models with —1.5 < n < 0, the simulations used 
were performed by S. White, with the same code as in 
EFWD. The total number of particles is 100^ and the 
force resolution equivalent to a Plummer softening param- 
eter e = L/2500, where L is the size of the computational 
box. The n — —2 run was performed by E. Bertschingcr us- 
ing an adaptive P^M code. It followed 128^ particles with 
e = L/2560. The CDM simulation denoted CDMl in Table 
1, is the same as in Gelb & Bertschinger (1994). It is a P^M 
simulation with 144^ particles, = 1, Ho = 50 km/s/Mpc, 
L = 100 Mpc and e = 65 kpc. It is normalized so that as (the 
linear rms mass fluctuation in a sphere of radius 16Mpc) is 
unity when the expansion factor o = 1. For CDM at early 
times, a < 0.5, we have used the results of a smaller box 
simulation performed by S. White. This simulation, denoted 
CDM2, has 100^ particles in a 50 Mpc box, with a force res- 
olution equivalent to e = 20 kpc. 

4.1 Numerical Resolution Effects 

In using N-body simulations to study gravitational cluster- 
ing, there are several numerical artifacts that need to bo 
carefully eliminated. In this sub-section, we focus on those 
effects that are most relevant in studying the small-scale 
dynamics of interest for this work. There are two critical 
factors: (i) The suppression of ^ on small scales due to force 
softening and 2-body relaxation, and, (ii) Effects duo to the 
initial conditions and the finite volume of the simulation 
box which limit the range of time outputs and length scales 
over which clustering follows the true dynamics. We provide 



estimates of the relevant numerical effects, and thus arrive 
at criteria for the range of scales on which the statistical 
measures used are reliable. For the scale free spectra these 
criteria are verified by testing for self-similar scaling. 

A comprehensive study and tests of numerical effects 
can be found in CBH, Baugh et al. (1995) and Gelb (1992). 
These authors have tested the effects of finite resolution, 
discreteness, and statistical fluctuations in the initial con- 
ditions on different statistical measures. Their results on 
small scales and on 2-point statistics arc relevant here, as 
the simulations used by them are of comparable resolution 
to those used in this study. Tormen et al. (1996) have made 
detailed tests of numerical resolution in the context of mea- 
suring density profiles of halos - these are also relevant as ^ 
is closely related to the halo profiles. 

Initial conditions. At early times the particle distribu- 
tion is affected by the cubical grid from which the particles 
are perturbed, and by transients due to the Zel'dovich ap- 
proximation used to set up the initial perturbations. The 
effect of the grid can be minimized by starting with "glass" 
initial conditions, as has been done in the simulations for 
— 1.5 < n < (see White 1994 for a discussion of "glass" ini- 
tial conditions). In general, to ensure that there are no arti- 
facts due to the initial conditions, the initial amplitude of the 
perturbations should be chosen small. A test for how small 
is sufficient is provided by checking for self-similarity, since 
the evolution of scale free initial spectra is expected to be 
self-similar once it follows the true dynamics. The initial am- 
plitude used in the scale free simulations with —1.5 < n < 
is such that the power on the Nyquist frequency of the parti- 
cle grid equals the white noise level (Efstathiou et al. 1988). 
For n — —2, as noted by Lacey & Cole (1994), departures 
from self- similarity persist for rather late times. Hence the 
initial amplitude needs to be lower for n = —2; it was chosen 
such that the dimensionless power on the Nyquist frequency 
of the particle grid is 4nk^P{k) = 0.03. 

Force softening. To suppress the effects of two-body re- 
laxation, force softening is implemented in computing short 
range forces in the particle-particle part of the siirmlations 
by using the equivalent of a Plummer force law: F{r) = 
Gm^r/{r^ + e^)^^^. The parameter e is given in Table 1 for 
the different simulations used. The results shown in Figure 
1 for (,{a,x) extend to 2e on the small scale end. As com- 
pared to the exact Newtonian force Gm? jr^ , the fractional 
error in F{r) at r = 2e is 28%. The results in Figure 1 show 
that for the first one or two points plotted, in the range 
2e < a; < 3 — 4e, ^ is underestimated, especially at ear- 
lier times. As discussed below, the cause for the artificial 
suppression of ^ at earlier times is usually due to 2-body 
relaxation and not force softening. 

Discreteness effects: 2-body relaxation. The effect of 

force softening is easy to quantify, as described above, and 
since e is constant in comoving coordinates, it is expected to 
affect the same range of scales at all times. But the distribu- 
tion at separations > 2e may still be inaccurate if randomly 
chosen particles do not have enough neighbours within such 
a separation, due to the effects of 2-body relaxation. A sim- 
ple estimate can be made because 1 -|-^(a;), when integrated 
over a sphere of radius x, is the mean number of neighbours 
within separation x from a randomly chosen particle. For ^ 
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Table 1. Parameters of the N-body simulations 



Simulation 


^ part 


Softening e/L 


PM Mesh 




n = 


100^ 


1/2500 


256^ 


2 - 


4 


n = -1 


100^ 


1/2500 


256^ 


2 - 


4 


n = -1.5 


100^ 


1/2500 


256^ 


2.5 - 


- 5 


n = -2 


128^ 


1/2560 


256^ - 432^ 


3 - 


5 


CDMl 


1443 


65Kpc/100Mpc 
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large enough this is essentially the total number of neigh- 
bours. Therefore with x = me ~ mL/2000, where m is of 
order 2-5, and N^^^-t/L^ ~ IQ^/L^ being the mean number 
density, we can estimate the mean number of neighbours to 
be: 

""•-5^(^)(t|.)(?)'- 

In obtaining the above equation, we have assumed that ^ cx 
cc~^ down to X = 0, with 7 = 1 — 1.8 for the spectra of 
interest. In the simulations however, ^ is suppressed on very 
small scales. The extreme case is given by setting 7 = 
which decreases Nnbrs by a factor of ~ 2. 

We use the empirical criterion that N„brs > 10 is re- 
quired to measure ^ accurately. The resulting minimum scale 
Xmin is giveu in Table 1 and is > 2e at the two earlier time 
outputs shown in Figure 1 (in stars and circles). It is typi- 
cally smaller than 2e for the last output, in which case we 
set Xmin = 2e due to the effect of force softening. Taking 
the minimum reliable scale to be the larger of the values 
that satisfy the criteria of accurate forces and sufficient par- 
ticle number successfully identifies the small scale regime in 
which the measured ^ departs from self-similarity. For CDM 
we have checked that the measured ^ in the two simulations, 
which have different values of Nnbrs at a given scale, are in 
very good agreement on the scales used in our analysis (see 
Figure 4). Our results for Xmin are similar to those of CBH 
after adjusting for the different particle number and force 
softening. 

Range of a;„;(a). At late times the nonlinear scale a;„;(a), 
defined by equation (p^, approaches the size of the box. 
This sets a maximum value of a beyond which the clustering 
is artificially suppressed due to the absence of modes with 
wavelength larger than the box-size. This problem is par- 
ticularly severe for spectra with n < —2 because they have 
more power on large scales relative to spectra with higher n, 
which in turns leads to a stronger nonlinear coupling with 
small scale modes. It is the principal reason for the difficulty 
in obtaining self-similar scaling for the n = —2 spectrum. 
For spectra with n > —2 we find that Xni{a) > L/10 is suf- 
ficient to obtain accurate self-similar scaling on the scales of 
interest. This is the criterion used to set the maximum time 
output for our analysis. 

At very early times there are transient effects from the 
initial conditions, and moreover clustering on small scales 
has not fully developed. We allow for an expansion factor 
of at least 4, and in the case of n = —2 of 12, between the 
initial time and the first output used for our analysis. At this 
time the nonlinear scale is in the range Xni — L/50 — L/IOO. 



Other effects. 

PM mesh: The mesh used to compute the long range forces 
in the PM (particle-mesh) part of the force computation 
also introduces an artificial scale. But if the particle-particle 
forces are computed for a conservatively chosen set of parti- 
cles at each timestep, as is done for the simulations used 
here, then the effects of the mesh are virtually absent 
(Bertschinger 1991). 

Poisson fluctuations: To compute ^ at a given radius, sums 
over particles within logarithmically spaced radial bins are 
made. Poisson fluctuations within each radial bin produces 
errors in the measured 5. This effect is also negligible as the 
number of pairs measured is large enough that the result- 
ing fluctuations in ^ are smaller than 1%, and would appear 
smaller than the symbols shown in Figures 1 and 2. 
Finite volume effects: At low wavenumbers the flnite number 
of modes present in a given bin in k causes statistical fluctu- 
ations in the initial power spectrum. As nonlinear evolution 
proceeds, these low— modes couple to the modes at higher 
k and can in principle lead to fluctuations in the distribu- 
tion at small scales as well. Also, modes with wavelength 
larger than the box-size are absent. In practice for the range 
of Xni{a) used, there is negligible effect on scales smaller 
than I// 10, which is the range of scales of interest here. The 
effect is not negligible on the higher moments of the den- 
sity, as emphasized by CBH and Baugh et al. But for the 
evolution of ^ the fluctuations between different realizations 
are about 2 — 3% on the small scales. We have checked this 
by comparing ^ for simulations started from three different 
realizations of the n = — 1 spectrum. 



5 RESULTS 

Self-similarity provides a powerful simpliflcation in study- 
ing the evolution of scale-free spectra in an Einstein-de Sit- 
ter universe. Using the self-similar solution for £^[a,x) from 
equation (^, we deflne a characteristic nonlinear scale Xni{a) 
as 

i{a,Xr,i) = l ■ (a) oca'/ (18) 

In the results presented in this section, the comoving dis- 
tance X will be scaled by x„;(a) at each time. Expressed 
as a function of x/xni{a) the correlation function ^, or any 
dimensionless statistic, is independent of time if the evo- 
lution is self-similar. This allows the measured correlation 
function at different times to be combined so that the effec- 
tive dynamic range of the simulation is extended. See Jain & 
Bertschinger (1996) for an analytical analysis, and CBH and 
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Figure 1. Sclf-similarly scaled correlation functions. Each panel shows ^{a,x) computed at four time outputs from simulations of scale 
free spectra with n = 0, —1, —1.5, —2 as indicated in the panels. It is plotted against x/x^^i (a), where ^(a, x„i) = 1 defines the nonlinear 
scale x„i{a) at each a(t). Between successive time outputs a;„((a) grows by a factor of ~ 1.5. The circles, stars, crosses and triangles 
show the outputs at successively later times, and therefore at later stages of nonlinearity. The dotted line is the linear ^. The dashed line 
shows the stable clustering slope for ^ ; its amplitude is given by the model of Sheth &c Jain (equations 19 and 20) . 
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Padmanabhan et al. (1995) for numerical tests of self-similar 
evolution. 

Figure 1 shows ^{x/xni{a)) for four different power law 
spectra. In each panel four output times are used to compute 
^. The figure shows that self-similar scaling is convincingly 
verified for all the spectra, with the possible exception of the 
n — —2 case which shows some scatter. Small departures 
from self-similarity can be seen at the smallest and largest 
scales. On small scales the departures arise because f is un- 
derestimated due to the effects of force softening and finite 
particle number. The effect is very small for all points ex- 
cept for the the very last, which corresponds to a separation 
~ 2e. On large scales the correlation function is affected, 
especially at late times, by the absence of waves that are 
larger than the size of the box. This causes ^ at late times 
to fall below the self-similar curve. 

As pointed out by CBH, a simple way to estimate the 
self-similar ^ for a given spectrum, is to take the largest 
of the different values measured at different times for each 
value of x/x„i{a). This should give a reasonable result for the 
cases where there is some scatter, as in n = —2, or the small 
scale end of all the spectra. The rationale is that numerical 
limitations cause a suppression of ^, both at small and large 
scales. However, we shall not need to do this because for the 
scales of interest, the measured ^ are sufficiently self-similar 
without any corrections. 

With the self-similarly scaled correlation functions for 
different spectra, the prediction of stable clustering can be 
tested. The dashed lines in each of the panels of Figure 
1 show the slope predicted by stable clustering. It is ev- 
ident that on the smallest scales, for about a factor of 5 
in x/xni{a), the measured slope of ^ agrees with the sta- 
ble clustering slope. As a first approximation this occurs for 
(, 3> 100 for all spectra. A more detailed examination of the 
slope of ^ in Figure 1 and of the mean pair velocity in Figure 
3 shows that the onset of stable clustering occurs at higher 
^ for spectra with larger n; thus for n = —1.5 stable clus- 
tering is valid for ^ ^ 200, whereas for n = 0, it is valid for 
^ > 500. 

Interestingly, the amplitudes of ^ which demarcate the 
onset of stable clustering correspond to approximately the 
same value of x/x„i independent of the initial spectrum: 
x/xni{a) ^ 0.07. This provides a convenient demarcation 
of the stable clustering regime as it is independent of the 
spectrum or of any feature of self-similarity, and can there- 
fore be applied to realistic spectra like the CDM spectrum. 
The variation with n in the value of ^ corresponding to 
x/xni{a) ^ 0.07 arises because spectra with larger n have 
steeper ^ in the regime of 0.07 < x/xni < 1. Since £,{xni) = 1 
by definition, ^(0.07a;„() increases with n. 

The amplitude of the stable clustering ^, i.e. the nor- 
malization of the dashed line in Figure 1, has been fixed 
by the analytical model of Sheth & Jain (1996). Their cal- 
culation assumed that on small scales ^ is determined by 
the density profiles of collapsed halos. They fixed the shape 
of the profile by assuming the validity of stable clustering, 
and used the Press-Schechter form for the number density 
of halos to then predict the amplitude of ^. The final result, 
equation (16) of Sheth & Jain, can be re-written as 

^{x/x„i) = C„ (—)"' , (19) 



where 7 = 3(3 + n)/(5 -I- n) is the stable clustering slope, 
and Cn is an n— dependent constant. For the spectra used 
in this paper it takes the values, 

C-a = {5.93,5.40,7.42,11.96} for 
n = {0,-1,-1.5,-2}, (20) 

respectively. 

The results show that again, except possibly for n = — 2, 
the measured ^ is in remarkably good agreement with the 
predicted value. For the range of small scales with ^ > 200, 
the measured ^ for —1.5 < n < is within ~ 10% of the pre- 
diction. For n — —2 the predicted value is larger by ~ 30% 
for the maximum measured 5. For this case, the large scat- 
ter between the different time outputs suggests that even the 
maximum measured ^ at a given x/x„i is an underestimate 
of the true value; the dashed line might therefore be accu- 
rate at the same level as the other cases. The verification of 
the amplitude of ^ predicted by equation (^9|) implies that 
both stable clustering and the Press-Schechter model pro- 
vide accurate descriptions of nonlinear clustering for a wide 
range of spectra. It also suggests that it is reasonable to as- 
sume that 5 is determined by the structure of dark halos on 
sufBciently small scales. The trend that the onset of stable 
clustering occurs for larger values of ^ as n increases is also 
in agreement with the estimate of Sheth & Jain, which as- 
sumed that only the half-mass radii of virialized halos have 
stabilized. Our results can be used to estimate the free pa- 
rameter in the spherical model for ^ used by Padmanabhan 
et al. (1995). The choice m = 3 in their equation (14) gives 
closest agreement with our results. 



5.1 Mean pair velocity 

The stable clustering hypothesis can also be tested by di- 
rectly computing the mean relative velocity of pairs of 
particle. In physical coordinates this velocity should ap- 
proach zero as the pair separation is made vanishingly small. 
In comoving coordinates therefore, the mean pair veloc- 
ity v{a,x) = ax must exactly cancel the Hubble velocity 
Hr = ax, so that the ratio —ax/Hr = 1. 

However v{a,x) is much more noisy than £^{a,x), be- 
cause of the high dispersion in the pair velocity. This is il- 
lustrated in Figure 2 where the measured —v/Hr is plotted 
against ^ for n — —1. It is only on averaging the data from 
three realizations that a meaningful result is obtained. At 
least five realizations for each spectrum would be needed to 
directly measure —v/Hr to better than 20%. In the absence 
of such extensive simulation data, we use an indirect method 
to measure —v/Hr. 

As described in Section 2, the evolution of ^ is connected 
with V via the pair conservation equation. Therefore to the 
extent that the simulation is accurate, v can be measured 
using the evolution of ^, and should give the same result as 
would be obtained from a simulation with enough particles 
to reduce the noise in the directly measured v. We use the 
form of the pair conservation equation given in equation (^ 
to solve for v. We evaluate the time derivative of ^ by tak- 
ing logarithmic derivatives. Since the evolution of ^ and ^ 
is close to a power law in a, taking logarithmic derivatives 
reduces the error introduced by the finite time interval be- 



Does Gravitational Clustering Stabilize On Small Scales? 9 





Figure 3. The ratio of the mean pair (peculiar) velocity to the Hubble velocity, —v/H7', as a function of the mean correlation function 
^. The pair conservation equation is used to solve for —v/Hr using the evolution of ^{a,x), as discussed in the text. Results are shown 
for three output times, corresponding approximately to the outputs shown in Figure 1, for the same set of scale free spectra. The near 
coincidence of the curves at different times verifies self-similarity, while the approach to the dotted line —v/Hr = 1 is consistent with 
stable clustering. 
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Figure 2. Mean pair velocity for n = —1. Simulations of three 
different realizations of the n = — 1 spectrum are used to directly 
measure —v/Hr. The average value is shown by the circular data 
points, with the error bars showing the l-cr deviations of the dis- 
tribution. The solid curve shows —v/Hr computed using the pair 
conservation equation as in Figure 3, for the same output time as 
the circular points. 



tween successive outputs. To estimate —v/Hr for a given x 
at output time Ui we use the finite difference equation 



log^(ai. 



log5(ai 



logai+i - loga-i-i 



(21) 



While the mean pair velocity contains the same informa- 
tion as the evolution of ^, it tests for stable clustering more 
directly for generic spectra like CDM which do not evolve 
self-similarly. 

The results for the inferred ~v/Hr are shown in Fig- 
ure 3 for three time outputs of the scale free spectra with 
n = 0, —1, —1.5, —2. Self-similarity requires that the data 
at different time outputs coincide, since both —v/Hr and 
^ are dimensionless functions of x/x„i{a). On scales with 
^ < 1 the linear theory relation —v/Hr cx ^ is verified. For 
larger values of ^, the ratio —v/Hr turns over after reaching 
a maximum value of ~ 2 at ^ ~ 10. For f > 100, the ratio 
approaches unity, though it appears to continue to decline 
with increasing ^ for a given time. For the range of ^ for 
which the slope in Figure 1 agrees with the stable clustering 
value, the measured —v/Hr = 1 to within about 20% for 
almost all the data. While this is a satisfactory consistency 
check, —v/Hr does not appear to flatten into an asymptotic 
regime up to the highest ^ measured. 

The results for —v/Hr on small scales are less reliable 
than the results for ^ because the former depends directly 
on ^{x), which is an integral of ^ over scales smaller than x, 
and is therefore suppressed even for x > Xmin- This artificial 
suppression largely accounts for the continued decline in the 



measured curves below the line —v/Hr — 1. We verified 
that if ^{x < e) is fixed at the stable clustering slope, then 
on scales in the range 2e < x < 4e, shown by the last 3 
points at each time output, —v/Hr would indeed be almost 
flat. Larger simulations, which accurately measure —v/Hr 
on scales significantly smaller than in our simulations, are 
needed to test its approach to the line —v/Hr = 1. 

While the evidence for the existence of an asymptotic 
regime in Figure 1 for ^ is still weak, the results for ^ and 
the mean pair velocity are certainly consistent with stable 
clustering on the smallest scales probed. These results are 
in agreement with those of CBH, who used tree-code sim- 
ulations in their analysis, with force resolution comparable 
to our simulations and < 1/4 the number of particles. Pad- 
manabhan et al. (1995) however concluded that for Q = 1, 
stable clustering is violated as they obtained —v/Hr > 1. 
A careful comparison of their results with ours reveals that 
the results from their P"^M simulations are consistent with 
ours, and with stable clustering. But the results for n = —2 
and n = — 1 from their PM simulations are inconsistent with 
ours at the maximum ^ measured in their simulations (typ- 
ically a factor of 2 — 3 smaller than ours). The reasons for 
this discrepancy are not clear, and merit a detailed compar- 
ison of clustering in high-density regions in PM and P^^M 
simulations. 



5.2 Evolution of the CDM spectrum 

The results for ^ and v for the standard CDM spectrum 
with (Tg = 1 at a = 1 are shown in Figures 4-8. Figure 4 
shows ^ measured at 4 different redshifts in the range 0.2 < 
a < 1. We use results from both simulations to maximize 
the range of scales on which ^ can be accurately measured. 
At early times, for a < 0.5, when the correlation length 
Xni < 5 Mpc, we measure ^ from the simulation with a 
smaller box-size, CDM2, as it has better resolution on small 
scales. For a > 0.5 CDMl is used as x^i becomes larger 
than L/10 for CDM2. We are thus able to measure ^ down 
to 2;„;/30 for a — 0.2, and a;„;/100 for a > 0.4. Since the 
results of Figure 1 show that x/xni < 0.07 is required for 
stable clustering, we have a sufficient range of x to test the 
CDM model in the range 0.2 < a < 1. 

Unlike the scale free spectra, for CDM it is not possi- 
ble to scale self-similarly and compare ^ at different times. 
However, if stable clustering is valid, then expressed as a 
function of the separation r = aa; in physical coordinates, 
£,{a,r) must grow in time as a"^ (see equation (^)) so that 
the mean number of neighbours oc a~^^(a,r-) remains con- 
stant. We therefore test for the constancy of a~^^(a,r) on 
sufficiently small scales in the two panels of Figure 5 for four 
values of the expansion factor in each panel. Since the onset 
of stable clustering should occur approximately at a fixed 
value of ^, the physical scale at which this occurs increases 
with a. Figure 6 shows the evolution of the correlation length 
Tniia) with a. The values of r = 0.07 r„j(a) for the four val- 
ues of a in each panel of Figure 5 are shown by the circles 
at the bottom. Thus for the curve extending to the smallest 
r, the onset of stable clustering is marked by the left-most 
circle. The values of r = 0.07 r„i (a) are 20, 180,700kpc for 
a = 0.2, 0.5, 1 respectively. The rapid growth of r„i(a) makes 
it difficult to assess for stability for non-power law spectra 
by comparing ^ at different times. But within the region of 
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Figure 4. Correlation function for the CDM spectrum. The solid 
curves show ^{a,x) computed from two simulations of the stan- 
dard Cold Dark Matter power spectrum at four output times, 
with a = 0.21,0.32,0.51, 1.0. The data shown in circles are from 
a simulation with a box size L = 50 Mpc, and those in stars 
are from a bigger box with L = 100 Mpc. To check for the ac- 
curacy of the measured the data at a = 0.2 from both sim- 
ulations are shown. The dotted line shows the observed galaxy 
auto-corrclation function with slope = —1.8. 

overlap to the left of the circles, the different curves are very 
close to each other. The results are therefore consistent with 
stable clustering, as verified also by the plot of —v/Hr in 
Figure 8. 

Both the amplitude and slope of ^ on small scales are 
very well approximated by the dashed line in Figure 5, which 
is the same as for the n = —2 spectrum in Figure 1. Thus 
its slope is —1 and its amplitude is set by equations (jl9|) 
and (12(1) with n = —2. For the lower panel in Figure 5, with 
a < 0.5, the measured slope of ^ is shallower, as it is sensitive 
to the part of the CDM spectrum with n < — 2. However the 
n = — 2 line is still an adequate approximation. We conclude 
that for CDM-like spectra the slope and amplitude of ^ on 
small scales is not very sensitive to the precise shape of the 
spectrum, and is approximated rather well by the stable 
clustering result for a power law spectrum with n ~ —2. 

Finally we test a parameterization of the evolution of 
^ commonly used in studies of clustering of galaxies and 
Lyman-alpha clouds at high redshift. The correlation func- 
tion is taken to have the form: 

^(a,r)=a^+^(^)"™ , (22) 

where 70 = 1.8 is the observed slope of the galaxy correla- 
tion function, and ro — 5/!.~^Mpc is the correlation length 
at a = 1. The parameter e = in the stable clustering 
regime and 0.8 in the linear regime. Our results show that 
there are two main reasons why such a parameterization is 
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Figure 5. Scaled correlation function for the CDM spectrum. 
The two panels show (^{a,x)/a^ for the CDM spectrum plotted 
against the physical length scale r = ax. The different curves 
should coincide in the stable clustering regime. The upper panel 
shows expansion factors a = 0.51,0.60,0.81,1.0, and the lower 
panel shows a = 0.21,0.28,0.34,0.43. The data for a > 0.5 and 
a < 0.5 are separated as the numerical resolution limit on small 
scales precludes their comparison over a common range in r in the 
stable clustering regime. The slope and amplitude of the dashed 
lines in both panels are the same as the stable clustering values 
for n = —2 shown in Figure 1. The x-coordinate of the circles 
marked towards the bottom is the scale 0.07r„((a); for given a 
stable clustering is expected for r < 0.07 r„; (a) as discussed in 
the text. 



inaccurate in the context of any CDM-like model. First, the 
growth of S,{a,r) in the intermediate regime is much faster 
that a^, so the stable clustering and linear regime do not 
give the upper and lower bound for the growth rate of 
Second, the onset of the stable clustering regime, and of the 
intermediate regime as well, occurs on a scale which grows 
rapidly with a. Therefore to parameterize the shape of ^ by 
one or two power laws which remain fixed in time is not 
accurate. For general use, following Hamilton et al. (1991) 
and Nityananda & Padmanabhan (1994), universal fitting 
formulae have been proposed for the nonlinear ^ which ap- 
ply to any smooth initial spectrum (see Jain, Mo & White 
1995 and Peacock & Dodds 1996). 
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Figure 6. Correlation length for the CDM spectrum. The 
measured correlation length in physical coordinates, defined by 
^(a, r„i) = 1, for the CDM spectrum is shown for 0.15 < a < 1. 
The scale 0.07r,ii(a), shown in Figure 5 to mark the onset of 
stable clustering, can be read off from this plot for a given a. 

However, if for reasons of convenience, a parameteriza- 
tion like that of equation (^2|) is used, the results of Fig- 
ures 5 and 7 show that e = is a good approximation for 
r/vni < 0.07, with r„i(a) shown in Figure 6. In the regime 
0.07 < r/r^i < 1 the growth rate varies rapidly, rising to 
e ~ 2 before approaching the linear theory value e = 0.8 at 
r > r„; . Thus in the intermediate regime between the stable 
clustering and linear regimes, e is typically larger than 1, 
and is significantly underestimated if it is taken to be 0. 

Figures 6 and 7 can be used to estimate e for a CDM- 
like model in which the growth of ^ is parameterized as: 
^{a,r) — a'^'^^ S^{r), with ^ being a function only of r in 
some narrow range of a. Figure 6 gives rni{a) for any desired 
a = 1/(1-1-2), and Figure 7 can then be used to read off e 
at the range of physical separations r of interest. For CDM- 
Uke spectra and with Q = 1, the results of Figure 7 are not 
sensitive to the shape of the spectrum and to the choice of a 
in the range 0.2 ^ a ^ 1. This discussion is of course subject 
to the uncertainties due to the possible biasing of galaxies 
relative to the dark matter. 



6 DISCUSSION 

The results presented in this paper can be summarized as 
follows. 

• The shape and evolution of the correlation function ^(a, x) 
is consistent with the stable clustering prediction of equation 
(|). We have verified this result with the power spectrum 
meeisured from the simulations presented here as well. 

• Direct measurement of the mean pair velocity is not suf- 
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Figure 7. Growth of ^(a, r) for CDM. For a power law parameter- 
ization of the form ({a, r) = a^^" ^(r), e is shown as a function 
of r/r„i. The circles and stars are for a = 0.28 and a = 0.51 
respectively. The variation of e with r/rni shows that a simple 
power law form for all r is not accurate, and that in the regime 
0.1 < r/r„i < 1, the growth rate is much faster than the com- 
monly used values of e = or 0.8. The value of e can be obtained 
as a function of r for any desired a = 1/(1-1-2) by using Figure 6 
to get r„i(a). 



ficiently accurate on small scales. We have therefore solved 
the pair conservation equation to estimate —v/Hr which ap- 
proaches unity on small scales, as required for consistency 
with the results for 5. 

• We find that the onset of stable clustering occurs at 
x/xni{a) = 0.07 for all spectra tested. This provides a use- 
ful way to demarcate the stable clustering regime for generic 
spectra. The range of ^ over which stable clustering is veri- 
fied is typically 200 < C < 2000; it is higher for initial spectra 
with more small scale power (n ~ 0). 

• For the CDM spectrum we find the range of scales for 0.2 < 
a < 1 for which the evolution of ^ is consistent with stable 
clustering. The combination of simulations with two differ- 
ent box-sizes enabled us to study clustering on very small 
scales, typically down to comoving scales = 40 — 120 kpc 
(with h — 0.5), while retaining the nonlinear influence of suf- 
ficiently long waves on the dynamics. The physical scale de- 
marcating the stable clustering regime rapidly increases with 
time, and is about r = 20, 180, 700 kpc for a = 0.2, 0.5, 1 re- 
spectively. 

• The amplitude of the nonlinear ^ agrees to better than 20% 
with the prediction of the Press-Schechter model of Sheth 
& Jain (1996) for nearly all the spectra studied. For the 
CDM spectrum the results agree with the predicted slope 
and amplitude of the n = —2 spectrum. The agreement of 
the model for different initial spectra is remarkable in view 
of the fact that their nonlinear behavior is quite different: 
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Figure 8. Mean pair velocity for the CDM spectrum. The mean 
pair velocity is computed using the pair conservation equation 
as in Figure 3. The three curves are for a = 0.3, 0.6, 0.8. Note 
that unlike Figure 3, for CDM-like spectra the curves at different 
times are not expected to coincide. They do however approach 
the stable clustering value —v/Ht = 1 for g > 200. 



e.g. for n = the nonlinear ^ is suppressed relative to the 
linear ^ by an order of magnitude, while for n ~ — 2 it is 
enhanced by about the same factor. 

• We emphasize that the commonly used parameterization 



,3 + e 



(r/ro) 



with e = is accurate only in the sta- 



ble clustering regime of r-/r„; < 0.07. It can severely un- 
derestimate the growth of ^ in the intermediate regime of 
0.07 < r/rni < 1 in which the growth of ^ varies rapidly 
with r, ranging from e ~ — 2. The results shown in Fig- 
ures 6 and 7 can be used to estimate e for given r and a for 
CDM-like spectra. 

For a numerical study aimed at testing an asymptotic 
regime it is appropriate to conclude on a note of caution. 
Small scale clustering is closely linked to the structure of 
the inner parts of dark halos. Both the number density and 
density profile of very small halos, with M <C M*/10, are 
poorly resolved in N-body simulations - this can crucially 
affect the shape of ^ in the small scale limit (S. White, pri- 
vate communication; Sheth & Jain 1996). 

Thus while the simulations used in this study provide 
strong evidence that gravitational clustering stabilizes on 
scales smaller than ~ 1/lOth the nonlinear scale, it is prema- 
ture to draw conclusions about a possible asymptotic regime. 
We have shown results typically over a factor of 5 in length 
scale, and an order of magnitude in ^, by using self-similar 
scaling to combine results from different output times. The 
competing demands of small scale resolution and a simula- 
tion box sufficiently bigger than the non-linear scale, do not 
leave enough dynamic range to see a convincing asymptotic 
regime on small scales (say a range of scales S> 10). 



We have carefully analyzed the effects of limited numer- 
ical resolution to rule out the possibility of a conspiracy of 
numerical artifacts over the range of scales that were used in 
our analysis. Therefore we can conclude with confidence that 
any departures from stable clustering could occur only on 
scales with ^ ^ 10*, which lie beyond our resolution limits. 
We shall have to wait for the next generation of high reso- 
lution simulations, with A'^ ^ 10* particles, to find possible 
new aspects of nonlinear clustering. 
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